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Fusion  of  Range  and  Luminance  Data 
&om  Laser  Radar  Systems 

Griff  L.BHbro 
November  8, 1993 


Abstraet 

In  thif  project,  a  technique  wat  developed  for  computing  an  op¬ 
timal  description  of  some  physical  scene  by  using  global  optimisation 
and  Bayesian  techniques  to  combine  correlated  outputs  of  diverse  non¬ 
linear  sensors.  Two  laser  range/intensity  imagers  were  considered,  but 
the  major  application  was  to  obtain  the  best  known  restorations  of 
very  noisy  PD-,  T1-,  and  T2-weighted  Magnetic  Resonance  Images  by 
fusing  the  information  in  all  the  input  images. 


1  Foreword 

Sensor  fusion  is  an  automatic  technique  for  describing  some  single  thing 
or  event  by  combining  relevant  data  &om  several  sensors.  Sensor  fusion  is 
becoming  more  important  as  more  sensors  are  integrated  with  computers. 
The  general  objectives  of  sensor  fusion  are  three: 

1.  To  reduce  noise, 

2.  To  increase  the  extent  or  detail  of  the  information,  or 

3.  To  tolerate  faulty  or  missing  data. 

This  report  is  restricted  to  an  imaging  application,  but  our  formulation 
was  consciously  developed  to  be  more  general  as  will  be  discussed  after  imag¬ 
ing. 
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1.1  Imaging  Application 

For  imapng,  the  objective  is  to  is  to  usefully  combine  several  nonlinearly  re¬ 
lated,  incomplete,  distorted,  or  noisy  signals  to  obtain  an  improved  restora¬ 
tion  or  reconstruction  of  the  true  scene.  This  imaging  problem  is  easiest  for 
pixel  registered  images,  which  in  this  case  is  multiple  image  restoration  to 
remove  noise  without  degrading  sharp  edges.  We  considered  three  imag^g 
applications  for  sensor  fusion. 

1.  The  Tri-Services  Laser  Radar,  a  time-of-flight  imaging  ranger  with  five 
channels:  coarse  range,  fine  range,  intensity,  doppler,  imd  passive  IR. 

2.  The  Odetics  range  camera,  a  triangulation  imaging  ranger  with  two 
channels:  range  and  luminance. 

3.  Commercial  msLgnetic  resonamce  imaging  system,  a  radio-frequency  im¬ 
ager  in  three  dimensions  controlled  by  user-specified  echo  time  TE  and 
repetition  time  TR,  with  three  channels:  proton  density,  single  nuclear 
relaxation  time  Tl,  and  nuclear  pair  relaxation  time  T2. 

We  formally  addressed  each  problem  in  the  first  year,  but  settled  on  the  MRI 
application,  item  3,  because  of  the  widespread  use  of  these  diagnostic  medical 
systems  and  the  resulting  availability  of  data. 

1.2  General  Formulation  of  Sensor  Fusion 

We  developed  a  general  formulation  of  low-level  sensor  fusion  as  minimization 
problem.  Find  the  optimal  estimate  of  an  unknown  measurement,  scene,  or 
dataset  /,  given  a  prior  expectation  V[f]  of  /  and  several  incomplete,  noisy, 
or  distorted  observered  images  Gm[f]  of  /.  We  model  /  as  a  2D  or  3D 
array  of  scalars  or  vectors.  Image  formation  for  channel  is  modeled  as  a 
functional  of  the  unknown  /.  V'e  considered  several  noise  models,  including 
Poisson,  Rayleigh,  and  additive  zero-mean  Gaussian,  although  the  following 
report  is  restricted  to  additive  zero-mean  Gaussian,  which  we  model  as  a 
random  image  for  each  channel. 

We  model  V  as  a  functional  of  the  unknown  /.  We  considered  a  smooth¬ 
ness  prior  in  terms  of  the  Lapladan  ^  ^  quadratic  vari¬ 

ation.  We  considered  a  flatness  prior  in  terms  of  the  sobel  I^.(Va/)^.  We 
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also  considered  priors  that  limit  the  search  to  a  spedfted  interval  of  values. 
The  following  report  restricts  itself  to  the  last  two. 

The  resulting  problem  is  then  to  minimise  the  objective 

=  E  -  c»(/)ll’  +  E 

m  p 


over  all  possible  values  of  the  measurement,  scene,  or  dataset  /.  The  fol¬ 
lowing  report  restricts  itself  to  the  problem  of  Magnetic  Resonance  Image 
Restoration. 


2  Statement  of  the  Problem  Studied 

We  studied  the  specific  problem  of  optimally  combining  several  nuclear  Mag¬ 
netic  Resonance  Images  (MR  images)  to  obtain  the  best  possible  description 
of  human  tissue.  We  formulated  this  multiple  image  optimization  problem 
mathematically  as  follows. 

Let  6  be  a  measured  vector-set  of  images 

G  =  {«c}c=l..  (1) 

where  d  is  the  number  of  channels  in  the  vector-set,  and  where  ge,i  represents 
the  c-th  channel  value  associated  with  the  t-th  pixel. 

Using  similar  notation,  let  S(F)  represent  the  undegraded  ideal  images  as 
a  deterministic  function  of  F  where  F  are  the  undegraded  ideal  basis  images, 
and  let  N  represent  additive  noise  such  that  G  =  S  -I-  N.  Note  that 

F  =  (2) 

where  p  is  the  number  of  basis  images  in  the  vector-set,  and  where 
represents  the  value  associated  with  the  *-th  pixd  of  the  V^'th  basis  image. 
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2.1  Bayesian  Model 

In  Bayesian  restoration,  the  most  acceptable  result  is  the  result  with  the 
highest  probability  of  occurrence.  Let  F  be  an  estimate  of  F.  Bayes’  rule 
gives  the  posterior  distribution  [8]  of  F  pven  the  data  G  as 

P(G1F)P(F) 

P(G) 

That  is,  the  conditional  probability  of  occurrence  of  a  specific  restoration  F 
pven  the  data  G  is  equal  to  the  conditional  probability  of  occurrence  of  the 
data  G  given  the  specific  restoration  F  times  the  probability  of  the  occurrence 
of  the  specific  restoration  F  divided  by  the  probability  of  the  occurrence  of 
the  data  G.  We  refer  to  P(G|F)  as  the  “noise  term”,  and  it  describes  the 
noise  distribution.  P(F)  is  called  the  “prior  term”  and  it  describes  the  a 
priori  distribution  which  can  be  chosmi  using  a  priori  knowledge  about  F. 
Obviously  P(G)  is  constant  and  independent  of  F,  so  in  order  to  maximue 
the  posterior  distribution,  we  need  only  maximise  P(G|F)P(F). 

2.2  Physical  Model  of  MR  Images 

The  function  S(F)  is  given  by  the  physical  model.  In  this  work,  one  simplified 
nonlinear  image  formation  model  [11]  is  used. 

=  Pi  exp  (-TEe/Tai)  {1  -  exp  (-TRe/r,<)}  (4) 

where  p,  Ta  and  Ti  are  basis  images  of  where  ^  =  1,2,3,  respectively. 
TEe  and  TR^  represent  the  echo  time  and  relaxation  time  used  during  acqui¬ 
sition  of  the  c-th  data  image.  Tx  and  Ta  are  nuclear  relaxation  times,  and  p 
represents  proton  density,  contributions  due  to  proton  flow,  and  MRI  system 
gain.  Most  brain  tissue  is  perfuse  with  slowly  moving  blood,  hence  the  data 
should  not  be  subject  to  large  variations  in  proton  flow.  This  work  does  not 
address  the  effect  of  proton  flow.  Our  data  was  acquired  with  MRI  system 
gain  held  constant  for  all  scans. 

Note  that  this  physical  model  is  undefined  and  exhibits  singularities  in 
the  gradient  when  Tx  or  Tj  equals  zero.  Tx  and  Tj  are  real,  positive  and 
bounded  below  in  time,  but  using  a  noninfinitesimal  step  size  during  gradient 
descent  requires  that  Tx  and  Tj  be  constrained  in  code,  otherwise  negative 


P(F|G) 
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values  of  T\  or  T2  might  occur,  causing  numerical  overflow.  Because  of  this, 
a  constrained  optimisation  technique  is  required  to  find  a  global  solution  in 
the  minimisation  process. 

2.3  Global  Optimisation  with  Mean  Field  Annealing 

Mean  Field  Annealing  (MFA)  is  based  on  Simulated  Annealing  (SA)  and 
derives  its  power  and  generality  from  that  popular  optimization  procedure. 
MFA  differs  from  SA  by  analytically  approximating  the  relevant  Gibbs  dis¬ 
tribution  rather  than  stochastically  simulating  it.  SA  works  by  gradually 
cooling  an  on-going  stochastic  simulation  of  a  Gibbs  distribution.  Mean  field 
theory  provides  a  deterministic  approximation  to  a  Gibbs  distribution  which 
also  can  be  cooled  in  the  same  way  to  produce  a  Mean  Field  Annealing 
(MFA)  algorithm.  Many  SA  algorithms  can  be  converted  to  analogous  MFA 
algorithms  that  run  in  1/50  the  time  required  by  the  SA  vernon[13,  1,  9,  2]. 
However  because  it  is  an  approximation,  MFA  does  not  inherit  any  guarantee 
of  convergence  even  when  the  analogous  SA  does  converge. 

In  earlier  work  we[l,  9,  2]  showed  that  simulated  annealing  could  be  ac¬ 
celerated  with  the  mean  field  approximation.  In  this  approach  the  important 
structure  of  P  is  approximated  with  a  more  convenient  distribution  Po  for  a 
sequence  of  falling  values  of  7.  In  this  section  we  provide  an  information- 
theoretic  procedure  for  studying  a  pven  difficult  P  using  an  essentially  arbi¬ 
trary  easy  Pq  by  minimising  the  entropy  of  Pq  relative  to  P,  or  equivalently, 
the  cross-entropy  or  Kullback-LeiblerflO]  distance  between  Pq  and  P.  This 
information-theoretic  procedure  leads  to  our  previously  successful  approach 
based  on  the  theoretical  tools  of  statistical  physics. 

Assume  we  have  another  positive  but  otherwise  arbitrary  distribution 
Po[s,  m].  It  is  useful  to  choose  Pq  to  be  easily  analyzed  and  to  have  adjustable 
parameters  represented  by  some  vector  m.  We  rewrite  Pq 

Po[s,m]  = -^exp(-l/o[s,m]/r),  (5) 

where  Zo  =  /dsexp(— 17o[s,m]/T)  which  in  general  depends  on  m  through 

Uo. 

The  entropy  of  Po  relative  to  P  is 

R  =  J  dsPo[s,m]  In  (6) 
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where  we  have  suppressed  the  dependence  of  R  on  the  vector  of  adjustable 
parameters  m.  Using  Equation  5,  we  rewrite  Equation  6  as 

R  =  j  -  In  Zo  +  U/T  +  In  Z).  (7) 

We  define  the  average  with  respect  to  Po  of  a  function  ^[s]  as 

(4,)  =  j  ds4,exp{-UolT)/Zo 

and  obtain 

«  = -i(Uo-U) -InZo  +  lnZ.  (8) 

We  define  Fq  =  — Tin  Zq  and  F  =  —T]nZ  and  obtain 

R  =  j{F<,-F  +  {U-  Uc)),  (9) 

It  is  known[10]  that  iZ[m]  >  0  with  equality  holding  if  and  only  if  Pq  =  P. 
Here  T  is  also  positive  so  that 

F<Fo  +  (cr-[ro},  (10) 

which  is  the  basis  of  our  mean  field  approximations  to  discrete,  continuous, 
and  even  problems  with  both  discrete  and  continuous  variables. 

The  mean  field  approximation  is  obtained  by  minimizing  Equation  9  with 
respect  to  m  to  find  the  tightest  bound  in  Equation  10;  mean  field  annealing 
involves  tracking  the  minimum  from  high  to  low  values  of  T[12].  In  the  case 
of  discrete  Sj  as  in  graph  coloring  or  binary  image  restoration[4]  it  is  useful 
to  choose 

Uo  =  -'^miSi,  (11) 

i 

but  in  the  present  context  of  problems  with  continuous  s,-  the  simplest  useful 
choice[9,  2]  is 

«>  =  5  P*<- ’»<)’■  (w) 

In  either  case  the  are  real. 
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3  Summary  of  the  Most  Important  Results 

We  have  discussed  three  important  results 

1.  The  formulation  of  low-level  sensor  fusion  as  an  optimization  problem 
that  admits  prior  expectation  in  a  rigorous  Bayesian  sense 

2.  The  application  of  Mean  Field  Annealing  (MFA)  to  the  problem  of 
sensor  fusion 

3.  The  development  of  a  specific  algorithm  with  produces  the  best  known 
restorations  of  Magnetic  Resonance  Images 


4  Publications  and  Technical  Reports 

This  project  has  resulted  in  four  publications. 

1.  A  journal  article  reporting  on  MFA  for  as  a  general  technique  and  for 
image  optimization  problems  speciiically[3]: 

Griff  L.  Bilbro,  Wesley  E.  Snyder,  Steven  J.  Gamier,  and  James  W. 
Gault.  Mean  field  annealing:  A  formalism  for  constructing  GNC-like 
algorithms.  IEEE  Transactions  on  Neural  Networks,  3(1),  1992. 

2.  An  conference  presentation  and  subsequent  proceedings  publication  of 
the  specific  MRI  problem  8tudied[7j: 

Steven  J.  Gamier,  Griff  L.  Bilbro,  James  W.  Gault,  Wesley  E.  Snyder, 
and  Y.  S.  Han.  M^netic  resonance  image  analysis.  In  SPIE  Proceed¬ 
ings  Vol  1904:  The  SPIE  and  1ST  Conference  on  Image  Modeling, 
1993. 

3.  An  journal  publication  that  formulates  the  specific  MRI  problem  as  a 
sensor  fusion  and  global  optimization  problem  and  reports  state-of-the- 
ait[6]: 

Steven  J.  Gamier,  Griff  L.  Bilbro,  James  W.  Gault,  and  Wesley  E.  Sny¬ 
der.  Magnetic  resonance  image  restoration.  Journal  of  MathemaUcal 
Imaging  and  Vision,  1993.  Accepted  for  publication. 

4.  A  second  journal  article  (currently  in  review)  discussing  refinements  of 
the  technique[5]: 
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Steven  J.  Garnier,  Griff  L.  Bilbro,  James  W.  Gault,  and  Wesley  E. 
Snyder.  The  effects  of  various  basis  image  priors  on  mr  image  map 
restoration.  Journal  of  Maihetnaiieal  Imaging  and  Vision,  1993.  In 
review. 

In  addition  we  are  preparing  a  clinical  article  for  the  medical  literature  in 
conjunction  with  radiologists  at  Bowman  Gray  School  of  Medicine. 


5  Participating  Scientific  Personnel 

Mr.  Stephen  F.  Gamier 
Mr.  Michael  McCormick 
Dr.  Griff  L.  Bilbro 
Dr.  James  W.  Gault 
Dr.  Wesley  E.  Snyder 


S.  J.  Gamier  earned  the  Masters  degree. 

S.  J.  Gamier  has  finished  his  PhD  thesis  research  and  expects  to  earn  his 
PhD  degree  early  in  1994. 
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